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?? ' Abstract 

A set of algorithms is presented for efficient numerical calculation 
of the time evolution of classical dynamical systems. Starting with 
^^ ■ a first approximation for solving the differential equations that has 

a "reversible" character, we show how to bootstrap easily to higher 
order accuracy. 
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1 The Problem 

We start by considering Newton's Law of motion in one dimension, written 
as a pair of first order time evolution equations. 

±(x(t)\( v(t) ) =M (x(t)\ (11) 

dt\v(t)J {f(x(t))J M \v(t)) 

where both x and v are time dependent functions to be determined at some 
later time t, given their values at some initial time t — 0. The force is given 
by some specified function f(x); and the quantity M is defined as the (non- 
linear) matrix/operator specified above. 
For simplicity I will write 

m = ( $> ) . (1.2) 

We assume that there exists an operator E(tM) that is the exact propagator: 

if>(t) = E(tM) V>(0), E(tM) = lim (1 + — M) N . (1.3) 

The addition property E(t 1 M)E(t 2 M) = E((ti + t 2 )M) follows. Alterna- 
tively, we may write, 

—E(tM) = M E(tM). (1.4) 

at 

This formalism is familiar in the case where M is a general linear oper- 
ator, and E is simply the ordinary exponential function; however, it is also 
appropriate for non-linear operators, as derived in reference [1]. 

Our objective is to show simple and accurate approximations to the op- 
erator E(5M), for small time-steps 5, for use in automated computations. 



2 The General Method 

Following the general method given in [1], we start by constructing an ap- 
proximate propagation operator R(S) with the following properties: 

R(5) R(-S) = 1, (2.1) 

and 

R(5) = E(5M + 5 3 X 3 + 5 5 X 5 + ...). (2.2) 



Rather than expanding the approximate result ip{t + 5) ~ R(S) ip(t) 
in a power series in 5, we represent R as the exact propagator for some 
other problem, which is expanded about the true one: based upon M. The 
restriction (2.1) means that only odd powers of 5 occur in the expansion 
(2.2). The quantities X3, X$, etc., are unknown. Our method will show how 
to eliminate those higher order errors, step by step. 

One more general property of the abstract propagator function E is the 
following. 

E(A)E(B) = E(A + B+ l -[A/B] + 1[(A - B)/[A/B]] . . .), (2.3) 

E(A)E(B)E(A) = E(2A + B--UA + B)/[A/B\\ . . .). (2.4) 

6 

This is the nonlinear extension of the Baker-Campbell-Hausdorff theorem 
for the product of exponentials of non-commuting linear operators. The 
only difference is that, instead of the commutator [A, B] = AB — BA for 
linear operators, we have the "slash commutator" [-A/.B] = A/B — B/A for 
nonlinear operators, as defined in reference [1]. 

Now we proceed. The initial operator R(5) is correct to order 5 2 and so 
we call it R 2 (S). Now we construct the following sandwich: 

R*(8) = R 2 (f35) R 2 ( 1 5) R 2 (f35); (2.5) 

and try to choose the constants /3, 7 so that 

R 4 (5) = E(5M + 5 5 Y 5 + ...). (2.6) 

By working with Equation (2.4) we find the simple requirements, 

2/3 + 7 = 1, 2/3 3 + 7 3 = 0, /3 = (2-2 1/3 )-\ 7 = -2 1/3 f3. (2.7) 

This new formula (2.5) may be read as follows: Take a step forward of 
length 1.351207. .. S, then take a step backwards of length 1.702414 ... S, 
then another step forward of length 1.351207. . .5. The result will be one 
step forward of length S - with errors of order 5 5 . 



3 Choosing R 

The real challenge now is to construct R(5), seemingly accurate only to first 
order in 5 but restricted by the requirement (2.1). 



Here is one suggestion, for the particular problem we started with (1.1) 
that is built in the "sandwich" manner. 



R 2 (S) = D x (5/2)D v (5)D x (5/2) 
x \ ( x + 8v \ 

v I \ V ) 




(3-1) 
(3.2) 



(3.3) 



\v + 8f(x) J' 

It should be apparent that this formulation is very easy to program for 
automated computation. On the other hand, it is rather cumbersome if 
one writes out explicit formulas for the overall result of this sequence of 
operations. 



4 Numerical Examples 

I have applied this method to a simple problem, the Kepler orbit in a plane. 
With the initial conditions x(0) = l,y(0) = 0, i; x (0) = 0,1^(0) = 1, I broke 
a complete orbit into N steps and saw what was the resulting error in y(N), 
which ought to return to zero. The results are shown in the tables below, for 
various values of N and various values of the source strength g (g=l gives a 
circular orbit). 

TABLES of computational errors 



Using R 2 


g=0.625 


g=1.0 


g=2.5 


N=100 


2XKT 1 


8xl(T 3 


2xl(T 2 


N=1,000 


2xl0" 3 


8xl(T 5 


3xl(T 4 


N=10,000 


2xl(T 5 


8xl0~ 7 


3x10^ 



Using i? 4 


g=0.625 


g=1.0 


g=2.5 


N=100 


3xl0- 2 


8xl0- 5 


2xl0- 3 


N=1,000 


3xl(T 6 


8xl0- 9 


2xl0- 7 


N=10,000 


3xl0^ iU 


8xl0" 13 


2xl0" n 



Each increase in the number of steps by a factor of 10 improves the 
accuracy by a factor of 10 2 if we use R 2 and by a factor of 10 4 if we use R4. 



Of course, R4 requires three times as many operations per step, compared 
to R2; but that seems a worthwhile investment since we can use many fewer 
steps for a given overall accuracy. 

One can readily extend this calculational method to involve any number 
of particles, with any forces. 

For comparison, I ran this same calculation using the popular Runge- 
Kutta method, at second order, and compared the results with those shown 
above for R 2 . Overall, one sees the same rate of improvement in accuracy 
as N is increased; and this is to be expected. For this particular problem 
I found that my method gave somewhat better accuracy at each level; but 
I would not offer that as a general rule without much further study; and 
I encourage others to try both methods on their own favorite problems. I 
will say, however, that the programming for my method was considerably 
simpler than that for the R-K method; and I expect that this aspect of the 
comparison is even more marked as one goes to the fourth order methods. 

What about the Richardson technique? As a general rule, if you calculate 
something with a small parameter 5 and know how it converges to the true 
answer as 5 — > 0, then you can accelerate convergence. For example, if you 
know 

A(5) = A + 5 2 X 2 + 5 4 X 4 + ..., (4.1) 

then you can do two calculations and combine the results as follows. 

\a{8/2) - l -A{5) = A + 5 4 Y 4 + .... (4.2) 

I used this method on the Kepler calculation, using R2, and found results 
that were slightly better than those obtained from using R4. This appears 
to be a nice alternative method. 

Next I added the Richardson extrapolation to the Runge-Kutta (second 
order) calculation of the same Kepler problem; and found that the results 
were far inferior to those just mentioned. There was some improvement 
in accuracy but significantly less than expected. The reason is that the 
asymptotic formula (4.1) is incorrect for the Runge-Kutta method: there is 
a term of order S 3 that belongs there. 

The main lesson from these experiments appears to be that the condition 
(2.1) on the lowest order approximation, which we might call "reversibility", 
is important. 



5 Velocity-dependent Force 

Here we start by considering a simple type of velocity-dependent force: 

d x 

-jp=f(x,v)=g(x)+vh(x). (5.1) 

Now, when we construct the approximate propagator R2($), the D x part will 
be the same as before but the D v part needs to be changed so as to guarantee 
that R(-S)R(8) = 1. 

The way we do this is to find the exact solution of the simple equation 

v = g + vh, (5.2) 

where we treat g = g(x) and h = h(x) as constants. The solution is easy: 

v (t)=v(0) + (v(0) + ^)(e ht -l) (5.3) 

ft 

and this shows us how to construct the operator D v . 

With this special result we can now address the case of a general f(x, v). 
We do this by taking another time-derivative of the original equation (5.1). 

x = v , v = w, w = g(x,v) + wh(x,v) (5.4) 

uf(x v) ofix v) 

9(x, v) = v ^ » h{x, v) = gj • (5.5) 

Now we construct the following. 

R 2 (5) = D x (5/2)D v (5/2)D w (S)D v (5/2)D x (5/2) (5.6) 

x \ ( x + Sv 
DA<))\ v \ = \ v | (5.7) 

w I \ w 



(5i 



x 

D „.(<)) \ v \ = \ r j . (5.9) 

w 

This technique also shows us how to handle the simplest first order 
equation x = f(x) by turning it into a second order equation, x — v and 







v = vf'(x), and then constructing R 2 as in (3.1). As another option, instead 
of using 

D V (S) :v^ve sr{x) (5.10) 

one could use 

«*> — £» 

6 Time-dependent Force 

Let's return to the original problem (1.1) and allow the force to be explicitly 
time dependent. 



(6.1) 



The natural guess is for the second order approximate propagator R 2 to 
be built from the original R 2 as follows. 

R 2 (S) = D t (S/2)R 2 (S)D t (S/2), D t (8) :t->t + 8. (6.2) 

This means that we proceed as before but evaluate the force /(£, x) with the 
variable t at the midpoint of each time interval. This formulation preserves 
the property (2.1). 

Then we can proceed to construct R4 just as before, using three of these 
operators R 2 with the same weight factors /?, 7 as in (2.5). 

7 Discussion 

The method described here appears to be a powerful, simple and versatile 
tool. 

It should be apparent how one can continue to improve the method, 
going from R4 to Rq , etc. While it is not easy to guess in advance what level 
of accuracy will be most efficient in any given problem, the programming 
procedures outlined above make it relatively easy to experiment and find the 
best approach. 



It is interesting that the equations (2.7) have another solution, one that 
goes into the complex plane, as follows. 

/3« 0.324 ±0.135i, 7 « 0.352 =f 0.270i. (7.1) 

Which is the best to use? My guess is that for conservative systems, the real 
solutions are best but for dissipative systems this complex scheme may be 
better. This whole area needs further study and experimentation. 
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